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The Kalman Filter is a widely used procedure in tracking algorithms. When 
normality assumptions are violated, the Kalman Filter performance tends to degrade. 
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I. INTRODUCTION 



Tracking algorithms are widely used in modern technology. Perhaps the 
most commonly used algorithm is the Kalman Filter (KF) (cf. [Ref. i| and [Ref 2] ) 
which is based on the following assumptions: 

a) The target movement model is known, and has normally distributed 

fluctuations. 

b) The sensors tracking the target provide unbiased, normally distributed 

measurements of current target "position”. 

If the assumptions are satisfied. KF is the best algorithm possible. In practice, 
however, the normality assumption may not be justified, and the KF performance is 
degraded. A number of modifications of the KF have been made by difierent authors 

based either on heuristics, analytic derivations, or both in order to improve KF 

performance in different situations (cf [Ref. 3] and [Ref. 41 ). 

In this thesis a tracking procedure is developed for the following one-dimentional 
basic (and classical) model: 



a) 

b) 



9 = 9 - 

,ir| n 



T 



S n-1 



T/t -r 1 J 1 1 J’ 



>"*); 



« = 0,1.2.... 20 

j — 1.2,3....m n = 0,\ 



SO 



where 

9, : -^ j is the target position after time step n 

y,_ , is the measurement of target “position ’ obtained bv sensor - at time steD 

-h 1 ./ 



In the model the target performs a random walk with normally distributed steps with 
mean zero and known variance a", and the / h sensor's measurement errors are 
unbiased and t-distributed with known parameters [ci ,$■) where a. is the degrees of 

J J J 

freedom and cr. is the scaie parameter of the t-distnbution. There are m sensors. 

This model violates the ciassicai KF assumptions by allowing thick-tailed 
outlier-prone) measurement errors. This model is designed for the situation where m is 
small, and it isn't possible to make use of the central limit theorem effect. 

The mathematical tractability of the model is complicated by the fact that there 
is no simple analytic form for an estimate of location of the target with Student-t 
measurement errors. 
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In Chapter II two possible ways to derive the classical KF equations are given. 
These derivations will help to motivate the procedure suggested in Chapter III. Some 
practical considerations regarding the implementation of the procedure will be given in 
Chapter III as well. 

In Chapter IV the simulation experiment implementing the procedure derived in 
Chapter III is described and numerical performance comparisons of the KF and the 
new procedure are presented. 

Chapter V contains final remarks about the new filter and further topics for 
future development. 

Appendix A contains tabulated results of simulation experiment. 

Appendix B contains details on simulation experiment implementation. 

Appendix C contains graphical presentation of the results from Appendix A. 
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II. KALMAN FILTER 



The basic model for the classical Kalman filter is that the target moves according 
to a random walk with normally distributed step sizes, l'hc measurement errors are 
also normal and statistically independent from target “position", i.e.; 



at 

b) 



0 . , = 0 - 

y _ . ■ = 0 
■'fiTlj rr 



°/I-r 1 ./’ 



<5 , . . '•'N(O.i’ ); 

/: + 1 .j ]" 



// = 0.1,2,. ..20 

j= [.2,2,...rn n — 0,1,2,. ..°0 



There are three alternative ways to obtain the Kl : equations: 

1. the maximum likelihood estimation approach. 

2. the minimum variance unbiased estimation approach. 

3. the minimization of the mean square of estimation error ■ 0 — 6 ). 

1 n iv 

Extension and synergism of the first two methods will lead us to the procedure 
described in Chapter III. 



A. MAXIMUM LIKELIHOOD 

A A 

Let 0 be the estimate of 0 after time step n and C be its variance. Since 0 is 

n n t n n 

normallv 1 distributed with mean 0 and variance C , the conditional likelihood of 

" n n 



eiven 0 and 



" i ^ n + hi n ■+• J t m 



ib; 



/\ > 

(6 - 
1 a + l n. 



(v , 

1 1 +■ l ,y n +■ 1 



L(0 ]9 ,y ,,y ) = e 



* *n 

7 = 1 



f.2.1) 



where v 0 = = C ■+■ f 

and v ; = <T ~ — Vanjy) (j = 1 m). 



'in the KF case normality is preserved at each estimation step as can be easily 
verified. 
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Finding the maximum for this likelihood is equivalent to finding minimum of the 
exponent, 
i.c. 



min 



I — — 

7=0 j 



(2.2) 



(dropping the first subscript on the j's and using y 0 = 0^ lor a compact notation). 
Expanding the squares we obtain 



2 



Q(9 ) = e 2 V --20 , 

^ n+t 'l + l y rt+1 U. ,17 

7=0 7 7=0 7 7=0 7 

fliis simple quadratic expression is minimized by letting 



( 2 . 3 ) 



(ff.l) 



a -+• l m 



<71 y 

V J. 
~o", 


9 n . ^,+ 1, 


+ 1 


m , 


, m , 


v I 


1 + V 


7 = 0 U j 





Since (AM) is a linear equation we may immediately obtain the variance of the 
estimate. 



V ar\ 



( $ z , v 2 



(K 2) 



A - 0 

C , = Vane ,) = 

n -r L a f i 



u . 






l 



If) (If 



7 = 0°; 



;=o y 7 



Vi 

7=0^7 



Equations (A. 1) and (A'.2) are the only ones needed to implement the KF. 

B. MINIMUM VARIANCE 

A 

Given y o =0 ({ and Vj~y nn . j ■ j= 1 m ; it is desired to get an unbiased 

minimum variance estimate (liMVE) of 0 , , . 

Since y' 0v yj,...^ are all unbiased estimates of 0, ; _ ( , any linear combination 



^ Vo - rU i J 'i + -- + a M 



m / m 



(2.4) 



constrained by ya.= 1 will also be an unbiased estimate of G 

“ 7 n 



+ 1 • 
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Since 



Var(n)=« 0 2 v’ () 4-« 1 2 v 1 + ...+a m 2 v m , 
finding the UMVIi is equivalent to solving the following problem: 



(2.5) 



min V = a ;) -r n + u 1 -v, + ... + «;v ffl 
S.T. a [r a r ...-a m =\ 



( 2 . 6 ) 



Given a Lagrange multiplier X, the Lagrangum ctjuation may be expressed as: 

~ = a o\ + V v l + + a J v m ~ + a i + - + 



(2.7) 



After taking partial derivatives and setting them equal to zero the resulting system 
of m ■*- 2 linear equations is: 

-x = o 



2 U , r, - a = 0 



2a v - X = 0 

m m 

- a n - a, - ... - a = - 1, 

0 1 m 



(2.S) 



Multiplying each equation bv the corresponding l/(2v.) and adding all the equations 
together results in: 



\ = 



V 1 

~ o u j 



Substitution tiien leads to 



(2.9) 



l 



a 

j 



' 1 
~0 U J 



( 2 . 10 ) 
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So, once again the final equations are equations (AM) and (A'.2) 



m y 

V J. 

A ~7 0 v. 

(K. 1 ) 0 = H = 

•'* + 1 m 

y - 



i 

(K. 2) C , = min V — 

at 1 m 

V I 

Note that the minimum variance approach does not require normality assumptions, 
but normality assumptions make it equivalent to the maximum likelihood approach. 
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III. W-FILTER 



Now to return to the basic model, target motion is a random walk with normal 
step size, but the m sensors have independent Student-t uistributed measurement 
errors. 

The basic model: 



a) 

b) 



y r. 



•1 = °n- 

1,/' 0 n + 



• 1’ 

d 



1/ 



c nTl ~N(0.x : ); « = 0 . 1 .:,...co 

1 7 ~ J= 1.2..1....//1 n =0,1.2.... 20 



It is assumed that our procedure is producing unbiased approximately normally 

A 

distributed estimates, i bis means that after time step //, 0 is approximately normally 
distributed with mean 0 and variance C . This is the same assumption as used bv 

n n k 

A 

[Ref. 3] and [Ref. 4], Based on this assumption we need to construct 0^ + ( and C + j . 

A. THE BEST 0 ^ . 

H 1 

1. The Criterion 

A 

Based on the assumption above, the conditional likelihood of 0 n+! given 0 f! 
a n d J,+ i,i--J’ N +i, m has thc form: 



A *■> 

<0 .-3 “ 

1 a ■+■ I t 



C r t" 



L (0 !0 



t - t* J n -t- l,r"”^T 1,; 



)= e 



•n 



c\d ) 

f 



d. + 1 

I 



/v — 0 1 

1 J n — l \ l 



( y *+\ 



(3.1) 



where c\ a) is some constant depending on the degrees-of-ircedom parameter d. of the 

j j 

y' 1 * sensor. Although it is possible to look for the maximum likelihood point, taking 
that approach has the following two interrelated potential problems. First there may 
be multiple local maxima present. In an extreme ease there may be as many as rn + 1 
local maxima. (See [Ref. 4j for local maxima conditions analysis when m= 1.) Second, 
in the multiple maxima situation, the global maxima may be tail and skinny, not 
having much likelihood to support it. One would not like to commit all belief on the 
basis of such evidence. 
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To deal with these problems use the following (smoothing) approach will be 
used: instead of choosing 0 n+l to maximize the probability of having 0 /l+1 in an 
infinitesimal probability element dO around 0 +1 , the estimator will be chosen 

so as to maximize the probability of having 0.,— j in an interval 

A A 

[0 n „ [ ~ w.0 n j_ j + m?J, the later interval having a finite (usually non-zero) predetermined 
length 2w. In other words the concern is not having an estimate exactly "on the spot” 
but rather having it within distance iv from true value of 0^ ( . 

In mathematical terms the approacii wiil be to find the solution to the 
equation: 



(IVU.a) 



max I(Q ,) = max 

n -t- 1 



3 . +• w 

T + 1 



• Q . - UJ 

'l -f 1 



uqe j , )di 

~ n n + 1,1 J n -t- 1,/n 



n -r l 



n -+■ 1 



In the equation (IV A. a) w is serving as a tuning parameter and has the 
following effect on the solution: 

When w= 0 solving the equation (If.l.a) is equivalent to finding maximum 
likelihood point. 

When w>0 and finite, equation (IV. l.a) tends to down-weight skinny peaks 
more than fat ones. 

When w is large enough, equation (If. l.a) will have a unique - solution. 
However it is generally neither required nor optimal to use very large 
values for w. 

A 

When w approaches infinity, any value off) + ( wiil satisfy the equation. 

The idea of using a non-zero w has additional appeal since, in some practical 
situations, occurrence of a large tracking error may degrade sensor performance for 
consecutive measurements. (Such sensor dependence on tracking performance is not 
reflected in the modei. ) 

2. The Solution Technique 

Finding a solution for equation (If. l.a) directly involves an exhaustive search 
with numerical integration; not a very exciting prospect at first sight. However it has 
the advantage of algorithmic simplicity and guaranteed global maxima. 



‘Since the likelihood function is positive and asymptotically approaches zero 
when 0 f;+ ( approaches ± oo this statement can be easily proved. 
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a. Analytic I iterative approach 

Alternatively, we may try to find an analytic solution to the equation 
( IV. 1. a ) by taking the derivative of /(Q^j) with respect to 0 /l+1 and setting it equal 
to zero. 



dm ,) 

9 t 



d0 



l 



[ rid- { 



' LU-0 >“ 
9 




cid ) 




1 



10 



.itl 



a 

e r 






e 



9 

C rf 



n. 




d; + l 



(3.2) 



Alter" setting the derivative equal to zero and performing canceilations 



(0 , + W -0 ) 



9 






9 

c 

9 




7 = 1 



L 






(y^u~ w ~ Q 



ifl 




i -t 



/\ a y 

(0 — i(J — 0 J 

1 9 l 9 







; = t 



l 




«./ 



A 

tu — 0 

9 >* l 



a 

7 



i - 1 



1 

d 



(3.3) 
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Next, taking logs of both sides of this equation results in: 



A a o 

2 C + c 2 



m d +1 



_ y 

j = i 



• in 



1+1 






(3.4) 



, (Q . — iv— s r 

1 ri •+• 1 n 

2 C + t 2 

a 



m d - \ 

y in 

7=1 " 



A 

v -r a> — 0 
1 ,; 



a 



n + l \ 1 

/ 7 



Multiplying through by ( — 1 expanding the squares, and pulling out the common 
denominator of the logarithm gives: 



A n <•) AO A A A ,\ 

0 2 + u;" + 0 2 + 29 w — 20 w - 29 ,9 m 

71+1 71 a+1 71 71 +It 



C +r 

a 



+ ^ id - l)/n(d o“ + fy — a> — 9 )“1 — Ind o 

— j 7 7 * *■ lj + 1 — 7 



7=1 



7 = 1 

(3.5) 



An 9 A o A A A A 

0 2 + u;" + 9 2 — 29 u; + 29 w - 29 9 gH 

T-t-l t /it 1 a a + l a 



C +r 

a 



+ y W — l)/n[cf o 2 + (y +u;— 9 ) 2 ]— V ; n <i 0 

— 7 7 7 * + lj « + 1 — 7 

7=1 ;= l 



Collecting terms and rearranging leads to 



49 71+ 1 U; 

C + c 2 



C 



o 

(f 0~ r- (v , — U) — 0 

9 \ ' , , w ; j ' + 1 

a r 1 Un : 

_ 1 1 a 

, = t (io'f (y 

J 7 7 n-lj 



19 ’i) "» 

_ V 



-717 — 9 t ,)“ 



(3.6) 



.And finally, 



\ -s c ~ z ~ m 

iWl.b) 9 , = 0 - -2 ^ 

1-<-l 1 4,y J 

7 = 1 



hi - 1 )ln 



d 0“ +<V . +177 — 9 ,)" 

J J ^ ri -r- i J 1 -I- l 

> A > 

a a“-r(v , -a>-0 J 4- 

j ; ^ a -t- 1 j a + 1 



Equation (IE. 1.6) is the first order condition for a local maximum for 
equation {IV. l.u). But perhaps a more natural way to check the first order conditions 

A ' A 

at point 0 /j + L is to see ii the likelihood function takes equal values at points 0 + t — vv 
and 0 



71 + 1 



u\ 
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When w approaches zero, equation (IVA.b) approaches the form: 



m (d.+l)(y .-9 ) 

A , _ 2 V 7 J + rc + l 



9 = 'e + (C +t“) 

n 4* l n n — ? 2 , / 

, = i a O -r ( V 

1 1 j j J 



A 

. -9 ,r 

1 -r lj a+l 



(3.7) 



(The right-most term is the first derivative with respect to 0,,4-j of the logarithmic 
function.) 

This is exactly the first-order condition for a maximum likelihood parameter value. 

Since 0 __ ( appears in the logarithmic term of equation (JIM. 6) 

is difficult to solve directly. One solution procedure is to use (ll'.l.b) in a recursive 

A 

manner to obtain 8 +l . The empirical experience with this approach indicates that 
unless is relatively large, direct application of equation (If'.l.b) does not converge* in 
most cases. 

b. Integral evaluation approach 

Instead of using a sophisticated iterative proceduie, the simple exhaustive 
search procedure described below will be used. 

Define 

A 

®min + 1,1 n + 1 ,ir) , o\ 

A A 

U =/AU7A*(U A’ 0-1 _U 1 ) 

max 1,/ir 

A 

The interval T=<9 . ,9 1 includes all reasonable candidate points for 9 _u . . It is 

possible to evaluate equation (IIM.Z?) at a finite number of points (say 20) from the 
interval T in accordance with the resolution required, and cake the one with smallest 

A 

violation of equation (JJ'.l.b) as the estimate 0.^ . . 

The approach of picking the point with the minimum violation of equation 
((T.l.£) docs not guarantee even locai maxima. However in some of the simulations 
that were performed it proved satisfactory unless there weie many cases when the 
likelihood function had large “Hat spots” and u> was small. This generally happens 
when the measurement variances are huge and/or the number of observers in is large. 

A procedure based on evaluation of first or higher order derivatives 
will be computationally costly and in general will be more difficult to solve because of 
local maxima and "flat spots” in the likelihood function. 



•'The "jumps" between iterations may be very large and the estimate may reach 
totally unreasonable values. 



Therefore, the safest procedure for finding the solution of equation 
{IV. La) is exhaustive search over the interval T with numerical integral evaluation. 
This is the procedure used in the simulation studies reported in the next Chapter. In 
the simulation, search and integral evaluation were performed in a single DO-LOOP 
(See Appendix 3 for details on implementation). The simulation was not written with 
computational efficiency in mind, nevertheless the computational burden did not reach 
an unacceptable levei. 

B. VARIANCE OF 0 , 

n ' [ ^ 

In order to keep the filter going once 0^, has been obtained, it is necessary to 
compute its variance C + j . Since equation (IK.l.u) is implicit and equation (IVA.b) 
is not linear, there is no natural analogy to equation ( K.2). 

I. Linearization approach 

/\ 

One possible approach to calculate the variance of0 <j+l is to approximate 
equation {IVA.b) with a linear equation obtained by a first order Taylor expansion of 

A 

the logarithmic terms around 0 >j+ j -y n J r t , for each j. Cancelling and rearranging 
terms results in the following approximation to equation (IVA.b): 

m 

(VV.l.c) 0 , =K~ l ( 9 + (C +c' 2 )V k v , ) 

n + 1 n n — a t 1 j 

7 = 1 

where 

d rl 

k - 

d o " -r w “ 

J ! 

m 

K = 1 + (C - c-’) V k 

n — j 

7 = 1 

This leads to: 

m 

C , =K~'\C ~ (C -Vf ^ k z v ) 13.11) 

n + 1 n. n. _ j j 

7=1 

where v . = Var(y> /i+ { J ). 



3 . 3 ) 



(3.10) 



21 



Empirical evidence indicates that in practice this approach performs very 
badly, giving values for C n which arc far too small. The C n depends only on n and 
converges to a value for large n. 

2. Minimization approach 
a. The Idea 

Having rejected the above linearization approach, a different one is 
introduced by the following motivation. 

.A A 

The "best'' estimate. , is some combination of m . 

In the case of the KF it is possible to pick a linear combination having minimal 
variance and remain consistent with maximum likelihood enteria. In the WF case 
this consistence no longer hoids, since eenerallv the estimate obtained bv finding the 
solution to the equation (IV. l.a) will differ from the minimum variance estimate. But in 
this case it is still possible find a minimum variance linear combination constrained 

A 

to vicld the "best” (alreadv known) value for 0 _ , . 

" " n i 



This means that m order to obtain an estimate of variance C ; ^_ j 

to solve the following constrained minimization problem: 

The symbolic formulation would be 

nun V = u ~i* 4- a, 4- ... 4- (/ - r 
0 0 11 m m 



S.T. 



A 

« n 0 4- a 

0 n 



'T/I T 1,1 



-a 



•r n ~ l 




+ a. -t- ... + a = 1 

0 I m 



we have 



(3.12) 



where 



L> 



o 




(3.13) 



u = Vany 

J n. + l 



i 




> = t.2 n 



(3.14) 



b. Formula derivation 

If m— 1. then no minimization is necessary and the problem has a unique 



A 

For notationai simplicity the first subscript will be dropped, and let j.v = 0 



solution. 



and 0 = O ;i+1 . Let Xj and X, be Lagrange multipliers. The Lagrangian formulation 
is then: 



se= 



m 

- V 



o 

cTu 
“ J J 
J = 0 






( 3 . 15 ) 



Taking partial derivatives* with respect to (L and X { , X ^ to obtain a system of m+j 
linear equations with m-i-3 unknowns gives: 



2v n a 



o^o 



2i\ a 



ni 



X jj*Q — 0 



2v a 

m rn 



’V m 




0 



V 

~o 



a v 

r j 



=3 



m 



-V = l 

7=0 J 



(3.16) 



Multiplying each one of the first m + 1 equations by the corresponding 1/v- and adding 
all together resuits in: 






m y 

V -J. 




7 = 0 



L’ 

J 



(3.17) 



Since ^a, = 1. 

Multiplying each one of the first m+ 1 equations by corresponding y.!v. , the (m-r2)" u 
equation by ( — 2) and adding all together produces: 



-\ 



9 

m y *• m 



U L “ 

7=0 7 J=0 



V 

" J 



V 

J 



A 
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(3.13) 
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After changing signs a system of two linear equations follows: 

Xj/1 + 20 

k l B+k 2 C=2 



where 



rn y ~ f Tl y 



j^o u j 



. v . 

7=0 J 



Applying Cramer s rule gives 



\ = 



2 ec - 20 

AC - fl 2 



c= 



- V 1 






(3.19) 



(3.20) 



Substitution vields 



\ = 



2 A - 20/3 
AC - B 2 



(3.21) 



a = 



Vr-s 



ZV 



(/ = 0,l,2,...,/n) 



(3.22) 



or exDlicitlv 



n r n y r n y m it 

y (gvI.VdL(Vi JYi 



1 - 



' = 0 y =0 / 



= 0 



= 0 



(/ . 



/ m V ^ i / m V . 

f V vi.( V 12 



7\ — y . — v 
7 = 0 7 



,=«V 7 



i3.23) 



Note that the denominator is equal to zero in ease — J.*j = ... = r . winch is a highly 
improbable event. 
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Since the I iessian matrix 



0 



2v n 






0 

0 



0 0 






r >?I 



is positive definite ( v.>0; for ail j), the solution is a global minimum. 
The final equation for tiie variance estimate is: 



<W. 2) C = 

n+ i 



m y~ m ^ 



V J. V _ _ V J. 

U J -o °J jTo U J 



m 

v 

jTo 



m m y 

, i h "V 1 - V a 



7»0 W ; 



. =0 U ) 



m y 

v ii 

~0 



y, 



- e 



* V 'J. 



j = 0 U j 



Now the equations (IV. La) and (IV. 2) contain all that is necessary to implement the 
WF procedure. 



C. 



1. 

■> 



THE WF PROCEDURE 



Find 9 . = miniv yy j. 

nun n^n — 1,1’ u n . 1 ,m 



0 = nia .\'( 0 .y 

max v n v . 

Fix the actual search resolution R (depending on T = |0 



i,r* 



•\?> 



1 ,m' 



• , 0 ]). 
nin max ' 7 



Set up the arid with resolution R over the interval 0 . — w, 0 +- w). 

1 . v mm ’ max ' 

(Let k = w: R) 



4. Evaluate the likelihood Function at each grid point. 

A 

5. Find the j by picking the "best” grid point solution of equation (IF. l.a) in 
interval T using sum of 2k -r 1 adjacent likelihood Function values For integral 
approximation. 

(See Appendix B for details on implementation). 

6. Compute C l by equation (fF.2). 



D. PRACTICAL CONSIDERATIONS 

A number of issues reiateu to implementing the WF procedure on a digital 
computer arc listed below. These considerations were used in in simulation described in 
the next Chapter. 

1. Grid setup 

The exhaustive search over the interval T = |0 .0 | mav be performed 

1 mm max' - 1 

over a pre-specilied number of points to be checked or using a pie-specified search 
resolution. 
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The first approach is wasteful when T is small and the second approach is 
wasteful when T is large. Thus, it would seem reasonable to specify both a desired 
resolution and an upper bound on the maximum number of points to be checked, 
which remain constant over time. 

2. Very dense observations 

In the case when T is verv short there mav be two complications: 

' A 

1) Because ot resolution problems 9 _ , mav take exact vaiue ot 9 . or 9 . 

In this case equation (IV.2) may force equal to the corresponding v. , 

wnich is not right. 

2) Computing C bv (fV.2) with v- values verv close together mav vield 
underflow and division by zero. 

In order to prevent these problems the following modification is suggested 
and was used in the simulations: if T is very short (say less than 3 times resolution ) 

A 

set 9 _ . to be eauai to midpoint of T and compute C , using (X. 2). 

3. "Very smaii” numbers 

When m is large, the likelihood function takes very low values. To avoid 
underflow, simply multiply it by a large constant. The constant that was used in the 
sunulation is 1024" 1 . 

One more modification was implemented to avoid arithmetic with very smaii 
numbers. If the exponent of the likelihood function first term was less than - 100. 
the likelihood function was set to be equal zero. The meaning of the modification 
is that the procedure considers any target step larger than 
'■'x (standard deviation) to be a practically an impossible step. (One must be careful 
with such modifications if large target lumps are possible, as couid occur if the model 
involved Student- 1 0-jumps). 
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IV. SIMULATION RESULTS 



A. SIMULATION DESCRIPTION 

A simulation experiment was performed on an IBM 3033 and 4381 at the Naval 
Postgraduate School. The programming language used was FORTRAN-77, and 
single precision was used. Normai and Chi-square random numbers were generated 
with the aid of IMSL library routines. Student-t random numbers were generated 
from normal and Chi-square random variables with appropriate degrees of freedom 
(see [Ref. 5] ). Tracking error statistics were computed using the HISTGP subroutine 
from the NON1MSL library at NPS. Figures were generated using an IBM 
experimental graphics package. GRAFSTAT. 

For each simulation replication the tracking sequence was performed for 100 
steps , (i.e. n+ 1 = 1,2,..., 100). The number of replications was 1000 in all cases. 

The KF and WF procedures (with different values for w) were carried out using 
the same random numbers. (With normally distributed target steps and Student-t 
distributed measurement errors). 

A 

Statistics were collected on tracking errors 0^^.^ for 

«+l = 1,25,50,75,100. 

The simulation was performed for 6 cases: 

Case 1 -one observer, < 7 , = 1 

1 

Case 2 -two observers, < 7 , = <7, = 1 

1 - 

Case 3 -three observers, n, = <7„ = <7.,= 1 
Case 4 -three observers. <7, = = 1: <7, = 3 

Case 5 -five observers. <* ^ = <7., = c> 3 = n_, = <7 C> = 1 
Case 6 -five observers, <7, =<7„ = <7,= 1; <7 . = <7- =3 

In all cases the measurement errors of the observers have the Student-t distribution 
with a’= 3 degrees of freedom and scaie <7., and the target has normally distributed 
step sizes, with standard deviation 7=1. In all cases the WF procedure was performed 
for values w= 0,1,2, 3. Excluding case 1, the desired grid resolution. DR, was chosen to 
be 0.1. In case 1, the desired resolution was equal to 0.05. Case 1 is most susceptible 
to problems resulting from a very short interval T described in the previous chapter. 
In ail cases the number of grid points on the interval T was limited by upper bound 
of 50. (See Appendix B for details on resolution and grid setup.) 



27 



After fixing a grid, the likelihood function was evaluated at each grid point over 
the interval of length T + 2 vv. The integral was approximated by the sum of 2k+l 
values of the likelihood function with k being the largest integer not exceeding vv/R. 

B. WF PERFORMANCE 

Figures C.l through C.6 display comparisons of the standard deviations 4 of 

A 

the tracking error. 9 . + , — 0 , . for the KF and the WF for cases l through 6, 

respectively. In ail the cases, the WF performance is considerably better than the KF 
performance. See Appendix A for tabulated resuits. 

Figures C.7 and C.S displav WF tracking error standard deviations with different 
values of w. The optimal w is somewhere in the vicinity of 1 or 2. but sensitivity to the 
exact value of w is low. See Appendix A for exact numbers. 

C. SENSITIVITY ANALYSIS 

For sensitivity analysis purposes, the tracking sequence was performed on each 
case with the same steps and measurements, but the WF procedure used the following 
(wrong) assumptions: 

ox = 0.5 or ox= 1.5 instead of <r.= 1.0 or or. = 3.0 (underestimatmg the 
measurement variance by a factor of 4) 

d. = 6 instead of T.= 3 (underestimating the measurement variance by a factor 
of 2) 

r = 0.5 instead of t= 1 (underestimating the target step size variance by a factor 
of 4) 

Additional runs were performed when target step size was a Student-t 
distributed variable having 3 degrees of freedom and variance equal to 1. 

In all sensitivity runs "he WF used w= 2 and the desired resolution was equal to 
0.1. Also, the maximal number of points in the interval T was 50. 

Figures C.9 through C. 14 display WF performance under sensitivity analysis 
conditions for cases 1 through 6, respectively. Clearly the WF procedure is very robust 



4 Ail the resuits presented in terms of sample standard deviations, i.e. a square 
roots of the unbiased estimates of population variances. Means are not presented 
since, in all cases for both filters, mean error values were very small. RMS (square 
Roots of Mean Square error) values were also computed in the simulation. In all runs 
inspected (the vast majority from ail runs) RMS values were were slightly smaller than 
the standard deviations values (due to division by N in RMS computation and division 
by N — I in standard deviations computation). 
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with respect to measurement error parameter estimation and less robust with respect 
to major alternations in target movement model. See Appendix A for tabulated 
results. 

The strange behavior of the standard deviation plot with t-distributed target 
step sizes is explained by the following phenomena: in very rare cases when the 

A 

estimate 0 and actual target position are extremely distant (the target makes a 

huge step), WF starts to “ignore' - ’ the incoming measurements, regarding them as 
outliers and sticks to the oid estimate. In such situation. C is increasing steadiiv but 
siowly (by It each time step - due to the wrong target model). After some time 
(which may be as many as 30 time steps), C becomes very large and WF stops 
disregarding observations and the estimate shifts toward the actual target position. 
(The last modification described in the previous Chapter was harmful in this 
situation.) 

Such undesired behavior may be treated by ad hoc adjustments to the WF 
procedure. For example, after detecting a steady increase in C. ; we may arbitrary set 
C n equal to a large number and recalculate the last estimate. This kind of adjustment 
may be proper when actually implementing WF if there is concern about the normality 
of target step size. None of these adjustments were implemented in the simulation 
experiment. It is anticipated that if an outlier-productive model of the target motion 
w r ere used, the effect noted would be automatically reduced, but to date no 
investigation has been made. 

Another comparison was made between the WF and KF procedures bv 
performing a simulation experiment with normally distributed measurement errors. 
The simulation was performed for cases la-da. Cases la-6a correspond to cases 1-6. 
respectively, with normally distributed measurement errors having the same 
variance as t-distributed errors in the original cases. Once again cne WF used w=2 
and the desired resolution was set to 0.1 with a maximum number of points in the 
interval T of 50. 

Figures C.15 and C. 16 display the comparison of WF and KF procedures 
for cases la-6a. It is evident that the advantage of KF over WF in cases la-6a is less 
than the superiority of WF over KF in cases 1-6. See Appendix A for tabulated results. 

In order to see how much measurement error assumptions may be violated, 
several more runs were made using the Cauchy distribution (Student-c with 1 d.f.) 
for generating measurement errors (WF assumed 3 d.f. and the correct second 
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parameter <J). In general the tracking was steady. However in very extreme cases, 
when the interval T was extremely long and the old estimate 0^ was not one of it's 
endpoints (actual resolution in such cases became extremely large), WF produced 
an unreasonable estimate and did not recover from the huge error. 

To overcome the problem, two alterations were made to the simulation 
implementation: 

The Maximum number of points to check was set to 500. 

In the case when ail grid points on the interval T corresponded to the integral 
having a value equai to 0. the new estimate was set to be equal to , the old 
estimate. The original lmoiementation oroduced 0 . for an estimate in such 

‘ min 

cases. (See Appendix 3 for a better approach by altering grid setup.) 

After making these two alterations, the runs were repeated. The tracking was steady 
with no special problems. The results are tabulated in Appendix A and may be 
considered very good. 

A selected number of cases was chosen to perform a limited variability 
demonstration. Simulation runs were made with 5 different pairs of seeds (one for 
target steps and one for measurement errors) for 1,2,3 and 5 identical observers and 
WF using vv=2. The resuits are displayed in Figure C. 17. Generally, variability 
between different pairs of seeds does not exceed the variability between separated 
time steps for a single seed. 
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V. FINAL DISCUSSION 



The WF procedure has a very general nature. It may be applied directly to 
any type of measurement error distribution. The distribution of measurement error 
must have finite variance and known density (even in tabular form). The WF filter 
performance with exotic distributions, using right or wrong distributional 
assumptions, needs to be investigated. 

(To appiy the WF procedure to models with non-normai target step size there 
must be an efficient way to evaiuate the conditional likelihood function, or else 

A 

additional assumptions have to be made regarding distribution of the estimate 0. t 
before and after movement of the target.) 

When applied to models with normally distributed measurement errors 
(and normally distributed target step sizes), the WF procedure (with any vv) is 
theoretically identical to the KF procedure. In practice, however, there might be some 
difference in performance due to finite resolution. 

WF has a built-in tuning parameter: the value of w. The optimal value for vv is 
definitely not zero, but the sensitivity of the procedure to the exact value of w is low. 
The best values of vv seem to be in the range (?. 2x), where t is the standard deviation 
of the target step size. 

Instead of using w as a static tuning parameter where it is kept constant, w may 
be computed dynamically, based on the latest observations. For example, w=T.A 
constrained to w>0.3 and w<J and the usuai procedure for grid setup (see Appendix 
B) produced very good results for all cases, but was not superior to the results 
obtained with constant vv. (Dynamic tv results were not represented in the previous 
Chapter). 

Another modification to be explored, if it can be justified by operational reasons, 
is use of a non-symmetric interval in equation •: IF. La), i.e. performing integration 
with lower limit 9 - w, and upper limit 9,. _ — vv-, . 

Equation (IV. 2) may be used not only as part of WF procedure, but as part of 

A 

any “good” procedure specifying 9 j and looking for C + 1 . 
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Optimal rule for dynamic w calculation may be discovered by further research. 
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The most noticable difference between the KF and the WF with respect to C n 
estimation is that in the KF case C n depends only on n, converging to some value as n 
increases, and in the WF case C n is more dynamic, depending on the actual 
measurement values. The dvnamic behavior of C reflects the variable amounts of 

* 71 

information received at each measurement step. This phenomenon gives an '‘inside - ’ 
performance measure on how weil the filter is doing. None of the cases checked by 
simulation required corrective actions due to abnormal C behavior (i.e. steady 
increase or large jumps), except when target step sizes were Student-t distributed. 

A iimited inspection of C values calculated by the WF procedure indicates that 

the C„ values computed are somewhat larger than the actual variances of tracking 

errors. It may conjectured that there is room for improvement in the WF 

procedure. On the other hand, forcing C n close to its true value (known from previous 

simulations) produces good but not superior results to those obtained in the regular 

way. using equation (IF. 2). This suggests that the WF procedure is fairly robust 

with respect to modest changes in the C n computation process. (A constant value for 

C mav be used if there is sufficient knowledge of the environment, i.e. the C values 
n ~ n 

that should be obtained are approximately known and no large jumps in C values are 
expected). 

Extending the WF procedure to the multidimensional case does not seem to 
pose conceptual problems. Equation (IV A .a) may be replaced by an integral over a k- 
dimensional rectangle. The analog to equation {W. 2) may be obtained by minimizing 
the determinant of the covariance matrix. However, the computational burden of 
solving equation {IVA.a) in the k-dimensional case and minimizing the determinant 
may be excessive. To actually implement the WF procedure in the 
multidimensional environment would require an efficient way to solve equation {IVA.a) 
and to minimize the determinant of the covariance matrix. 
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APPENDIX A 
TABULATED RESULTS 



This Appendix contains tables of tracking error standard deviations. The 
tables are organized by observers configuration. 

1. Table 1: One observer. 

2. Table 2: Two identical observers. 

3. Table 3: Three identical observers. 

4. Table 4; Three observers, one of whom is 3 times less accurate than the others. 

5. Table 5: Five identical observers. 

6. Table 6: Five observers, two of whom is 3 times less accurate than the others. 
The following notation is used for specifying target step size and measurement 

error distributions: 

• N(v) - normal with mean 0 and variance v 

• t(d,s) - Student-t with d degrees of freedom and scale parameter <7 = s 

When two types of observers are involved, the parameters are given for the more 
accurate ones. 

All cases above the double line use precisely the same target steps and 
measurement errors: only the filters are different. 

Filter notations: 



- KF 

- WF iv = i 
WFvv= flics’ 

- Wfvv= \,2d 



Kalman filter using correct variance 

W-filter using 3 d.f., -v = i and correct variance 

W-filter using 3 d.f.. sv=i and assuming the scaie of the t is 'icr 
rather than the true parameter <7 for each observer 

W-filter using w-i and assuming that the Student-t degrees of 
freedom are 2 d rather than the true value of d 



WFw= l.i it W-ffiter the correct measurement parameters and assuming the 
standard deviation of the target step size is 'it rather than the 
true vaiue of t 
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TABLE 1 
ONE OBSERVER 



Mvmnt 

step 


Measur. 

Error 


FILTER 


Track Length 

1 25 50 75 100 


N(_l) 


tC 3,1) 


KF 

WF vv = 0 
WF w= 1 
WF ;v = 2 
| WF iv= 3 


0.830 
0.842 
0.832 
1 O.S35 
0.837 


1.108 
1.1 o0 
1 1.096 
1.024 
| 1.033 


1.079 
1.082 
1 .037 
0.998 
1.021 


1.121 

1.176 

1.124 

1.040 

1.060 


1.117 
1.094 
1.048 
0.993 
1.031 i 




- - 


i WFvw=2.‘-i<r 
I WFw= 2.2 d 
WFiv= 2. !; r 


0.846 

0.821 

0.915 


1.017 

1.013 

1.146 


1.003 
1 0.993 
| 1.192 


1.042 
1 1.038 
1.209 


0.993 i 
0.998 
1 1.221 


t( 3..58) 
N(D 


1 WF w = 2 

t(l-I) WF w— 2 

' 


0.853 

0.8S4 


1.632 i i .702 1.141 

1.349 1.445 ; 1.362 


1.125 ! 
• 1.406 j 


N(l) 


N(3) KF 

WF w= 2 


0.857 

0.871 


1.203 , 1.166 | 1.124 
1.250 1.209 i 1.185 


1 1.131 ; 
1.192 



TABLE 2 

TWO OBSERVERS 



Mvmnt 

step 


Measur. 

Error 


FILTER 


1 


25 


Track Length 
50 ** "5 


100 




N(l) 


*43.1) 


KF 


0.735 , 


0.891 


i 0.859 


0.888 


0.879 






WF w = 0 


0.672 : 


0.S01 


1 0.757 


0.784 


0.771 


" 




WF v = I 


0.670 


0.792 


0.7*46 


0.7"*** 


0. 66 






WF w= 2 


0.682 


0.787 


0.743 


0.77*7 


0.769 : 






WF vj— 3 


0.713 


0.820 


0.779 


0.816 


0.808 i 




_ 


WF;v=2.-,tf 


0.084 


0.786 


0.7*46 


0.776 | 


j 0.767 






WFvv = 2.2 d 


0.6S3 - 


0.796 


. 0.749 


0.784 


0. "*5 


— 




WF'.V = 2.;-2 c 


1 0.835 i 


0.909 


: 0.902 


0.933 


0.917 1 


l ! 


t(3..5S) 




WF iv — *> 


1 0.7-9 


1 


1.499 


0.838 


1 0.857 


Nil) 


t(l.i) 


WF vv= 2 


\ 0.813 


LO04 


1.080 


: 1.006 


1.023 | 


N{1 ) 


N(3) 


KF 


■ 0.750 


0.924 


. 0.955 


| 0.908 


u.914 






WF iv = 2 


0.778 


0.959 


0.989 


0.946 


0.957 
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TABLE 3 

THREE IDENTICAL OBSERVERS 



Mvmnt 

step 


Measur. 

Error 


FILTER 


Track Length 

1 25 50 ~ 75 100 


Nfl) 


t(3,l) 


KF 

WF w=0 
WF w= 1 
WF w = 2 
WF 3 


0.663 

0.587 

0.586 

0.600 

0.040 


0.774 

0.682 

0.673 

0.669 

0.694 


0.748 

0.660 

0.653 

0.654 

0.687 


0.746 

0.633 

0.632 

0.645 

0 . 6 S 0 


0.778 j 
0.667 j 
0.658 
0.662 
0.707 I 




- - 


i WFvv= 2, ! i (7 
! WFw= 2\2d 
j WF)V= 2 ! 1 -b T 


0.602 

0.601 

0.778 


0.668 

0.674 

0.775 


i 0.658 
0.659 
0.802 


0.646 

0.652 

0.793 

i 


0.662 ! 
0.673 
0.820 


t(3,.5S) 

N(l) 


t(U) 


WF w = 2 
WF w- 2 


0.661 

0.750 


1.073 

0.843 


1.411 

0.890 


0.720 

0.859 


0.709 

0.836 


N(l) 


N(3) 


I KF 
WF w=2 


0.723 

0.745 


0.815 

0.849 


0.7S9 

0.809 


0.794 

0.826 

| 


0.768 

0.796 



TABLE 4 

THREE OBSERVERS, ONE 3 TIMES WORSE 



Mvmnt | 
step 


Measur. 

Error 


FILTER 


Track Length 

1 25 50 ~ 75 100 


NYl) 


t(3,l) 

X 


KF 

WF w = 0 
WF w= 1 
WF w= 2 
WF w= 3 


0.721 

0.642 

0.641 

0.659 

0.696 


| 0.887 
1 0.793 
1 0.778 
0.768 
0.795 


0.866 

0.766 

0.758 

0.736 

0.782 


0.871 

0.762 

0.754 

0.716 

0.784 


0.874 

0.769 

0.760 

0.729 

0.802 






WFw= 2.h<r 
WlTu= 2 2d 
WF;v=2.!h 


O .606 

0.659 

0.825 


0.7% 

0.775 

0.878 


0.774 

0.75b 

0.905 


0.778 

0.762 

0.897 


I 0.773 
0.768 ! 
0.923 | 


t(3..5S) 
1 N(l) 

i 


t(l".l) 


I WF vu — 2 
WF w = 2 

| 


0.740 | 1.079 
0.807 | 0.963 


j 1.455 
! 0.989 


0.809 i 0.805 ; 
0.937 | 0.958 ’ 


NT 1 ) 

[ 


N(3) 


KF 

WF ;i/= 2 


0.770 1 0.914 I 0.908 
0.798 ! 0.958 1 0.937 


0.900 

0.948 


■ 0.878 
j 0.914 
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TABLE 5 

FIVE IDENTICAL OBSERVERS 



Mvmnt 1 Measur. 

step 1 Error 

* 


FILTER 


Track Length 

1 25 50 75 100 


N(J) 


t(3,l) 


rr t" 

Kr 

WF iv =0 
WF vv = 1 
WF vv = 2 
WF vv =3 


O.060 
0.476 
0.480 j 
0.496 
0.547 j 


0.656 

0.536 

0.o35 

0.o37 

0.577 


0.616 

0.52S 

0.524 

0.530 

0.667 


0.632 

0.534 

0.532 

0.545 

0.588 


0.603 

0.512 

O.507 

0.513 

0.588 






WFvv= l.WG 
WFiv= 2.2 d 
WFvv= 2.1 = c 


0.508 

0.497 

0.681 


0.555 

0.542 

0.629 


0.546 

0.53^ 

0.640 


0.560 

0.552 

0.677 


0.535 

0.517 

0.631 


t(3..58) 

N(l) 


t(U) 


WF ;v= 2 
WF w= 2 


0.524 

0.671 


1.520 

0.668 


1.396 

0.689 


0.552 

0.683 


0.527 

0.655 


1 


N(3) 


KF 

WF w= 2 


| 0.618 
| 0.635 


0.657 

0.695 


0.663 
| 0.686 


0.650 

0.681 


0.655 

0.675 



TABLE 6 

FIVE OBSERVERS, TWO 3 TIMES WORSE 



Mvmnt 

step 


Measur. 

Error 


FILTER 


Track Length 

1 25 50 75 100 


N(l) j 


t(3,l) 

- - 

1 


KF 

WF vv = 0 
WF w= 1 
WF w = 2 
WF w =3 


0.645 

0.57 

0.567 

0.586 

0.617 


I'M 1»0— < 

oo-r-r-n- 

t'-vOKOvoO 

OOOOO 


0.739 

0.641 

0.633 

0.639 

0.670 


0.751 

0.647 

0.645 

0.654 

0.697 


0.700 
1 0.599 
0.601 
0.611 
0.649 


1 

| 


WFiv = 2.!:><T 0.606 

WFvv= 2.24 i 0.584 
WFvv= 2. '.2: 0.760 


pop 

b\ 

L- i -. O 


0.657 I 
; 0.644 
0.761 

; ... . . 1 


0.677 ! 0.636 ; 
0.664 1 0.616 
0.793 | 0.749 

1 1 


t(3..5S) : Vv'F w=2 

N(l) 1 td.l) j WF vv= 2 


0.635 

1 0.802 


; o.~20 
0.80 1 


i 1.435 
0.830 


0.690 1 0.630 | 
1 0.828 | 0.852 j 


Nil) M3) KF 

1 WF vv= 2 

1 


0.694 
s 0.724 

1 


■ 0.7" 1 
0.324 


0.790 

0.817 


0.782 0.772 : 

0.827 : 0.803 

1 1 
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APPENDIX B 

DETAILS ON IMPLEMENTATION 



In this Appendix the implementation of the simulation experiment is discussed. 
(There are many different ways to implement the WF procedure.) In particular the 
implementation of the grid setup and search for maximal integral over interval with 
length 2vv will be explained. 

I. GRID SETUP 

Two parameters (excluding w) control the grid setup: 

1) DR - the desired resolution (usually DR=0.1) 

2) NMAX - the maximal number of points to search ( usually NMAX = 50) 

The first step is to find interval T = [0 • .9 m 1, where 

1 L nrnrr. max J? 

®min = l, l »"—>■« -I- 1 ..>*)’ 

After finding T, the desired number of search points, n is computed: n = (T/DR)-f- 1. 
If n> NMAX, then n is set to be equal to NMAX. 

In the next step, the actual resolution R is computed: R = T/n . 

After fixing R. the number of points on interval w k. is computed: k=[w/R|. (If w is 
computed dynamically its value may change). 

Now the grid is fixed to have total numoer of n-f2k+I equally (R) spaced points 
having n+ 1 points on interval T itself (one on each end). 

The likelihood function values are evaluated for each grid point and stored in an 
array. (In the same order that the grid points appear on the real number line). 

A more economical way to set up the grid would be as follows: 

First, evaluate the likelihood function at the endpoints of interval T. If the 
likelihood is positive at both ends, proceed as described above. If the likelihood is 
essentially equal 0 on one or two endpomts, then eliminate the corresponding 
observation, recompute interval T, and repeat the procedure. (This was not 
implemented in the simulation experiment.) 
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2. SEARCH FOR MAXIMAL INTEGRAL 



Once all likelihood values are stored in the array, the search for the maximal 
integral was implemented in the following way. 

Initiation step: compute sum of first 2k -M values from the array. This sum is 
the integral approximation corresponding to 9 . . This is the first candidate for 




Search loop: advance on the array each time adding a new element to the sum 
and subtracting the oidest one. If the new sum is larger than the largest among the 
oid ones, take the corresponding grid point as the best candidate so far and record the 



largest sum. 

After n-r 1 repetitions the procedure gives the best candidate among the n— 1 grid 
points on interval T (and the corresponding integral value). 

This implementation approach was chosen mainly because of its algorithmic 
simplicity. Adding and subtracting REAL numbers repetitively might introduce 
some round-otf error. This is not a real problem in present case because: 

• 50 additions and subtractions is not a very large number. 

• The interest is in the grid point, not the corresponding integral value. 

• In close cases when round-otf error is reversing, the test can't be sure anyway 
because of finite resolution. 
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APPENDIX C 

COLLECTION OF FIGURES 



In this Appendix seventeen Figures presenting the WF performance (compared 
with the KF performance) are collected. 

« Figures C.l through C.6 present WF vs KF comparison for cases 1 through 6 
respectively. WF performance is given for values of w=0, 1.2,3. The cases 
correspond to following observer configurations: 

One observer 

Two identical observers 

Three identical observers 

Three observers, one of whom 3 times less accurate than the others 
Five identical observers 

Five observers, two of whom 3 times less accurate than the others 

• Figures C.7 and C.8 present WF performance with with different values of vv 

Figure C.7 - cases 1,2,3 
Figure C.8 - cases 4,5.6 

• Figures C.9 through C. 14 present sensitivity of WF performance with respect to 
model violations. Figures C.9 through C. 14 correspond to case 1 through case 6 
respectively. 

• Figures C.15 and C.16 present comparison of WF and KF performances for 
normally distributed measurement errors. Comparison presented for cases la-6a. 
Cases la-6a correspond to cases 1-6. but have normally distributed 
measurement errors with the same variance as Student-t distributed errors in 
original cases. 

Figure C.15 presents cases la-3a 
Figure C.16 presents cases 4a-6a 

• Figure C.17 displays variability of the simulation using different seeds. Five 
pairs of seeds are used. Variability is presented for cases 1,2, 3. 5. 



39 



0.8 



CM 




40 





41 




42 




0.2 0.4 0.6 0.8 




43 



0.2 0.4 0.6 




44 




Figure C.6 WF vs KF comparison for Case 6 - Five diderent observers. 
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Figure C.7 Sensitivity to w for Cases 1-3. 
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Figure C.8 Sensitivity to w for Cases 4-6. 
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Figure C.12 WF Sensitivity Case 4 - Three Different observers. 
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Figure C.16 WF vs KF comparison for Case 4a-6a - Normal measurement errors. 



55 



o 



oq 

o 



CD 



O 



CN 

o 



O 





o Case 1 

* Case 2 

x Case 
Case 



TRACK LENGTH 

i i l l i i 1 I l I 1 

0 20 40 60 80 100 



Figure C. 17 Variability of results. 
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